#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 23 01:48:35 2018

@author: Alexandre Dauphin

This script plots the Figure S6 of the supplementary material of the paper.
"""

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import brentq
from matplotlib.colors import to_hex


#### Definintion of the colors of the curves ####

n=np.array([20,50,100,200,400,800])
tn=np.size(n)

colors = plt.cm.winter(np.linspace(0,1,tn))
colors=colors[::-1]


#### Boundaries in the thermodynamic limit of the winding number, obtained ###
#### by the divergence of the localization length, as found in Ref. [5]. ####
m=1.12

def l(w,m):
    w1,w2=0.01,w
    f=np.abs(2+w1)**(1/w1+0.5)*np.abs(2*m-w2)**(m/w2-0.5)
    g=np.abs(2-w1)**(1/w1-0.5)*np.abs(2*m+w2)**(m/w2+0.5)
    z=np.log(f/g)
    return z

vecd=np.arange(0,6,0.001)
td=np.size(vecd)
wmin=brentq(l,0.1,3,args=(m))
wmax=brentq(l,3.,6,args=(m))
wind_inf=np.zeros(td)
wind_inf[wmin<vecd]=1
wind_inf[wmax<vecd]=0


#### Plot of the winding number for different system sizes ####

plt.figure()
plt.plot(vecd,wind_inf,'--',color='#989898')
plt.plot([0,6],[0.5,0.5], color='#d3d3d3')
for i in np.arange(tn):
    a=np.load(str(n[i])+'.npz')
    vecdis=a['vecdis']
    matwind=a['matwind']
    plt.plot(vecdis,matwind,color=to_hex(colors[i]),label=str(n[i]))

plt.xlabel('W')
plt.ylabel('Winding '+r'$\nu$')
plt.legend(loc=2)
plt.xlim(0,6)
plt.title('Figure S7')

plt.savefig('figure_s7.pdf')